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Abstract 

We apply the density of states approach to the Z 3 spin model with a chemical potential p. For determining 
the density of states we use restricted Monte Carlo simulations on small intervals of the variable for the density. 
In each interval we probe the response of the system to the variation of a free parameter in the Boltzmann 
factor. This response is a known function which we fit to the Monte Carlo data and the parameters of the 
density are obtained from that fit (functional fit approch; FFA). We evaluate observables related to the particle 
number and the particle number susceptibility, as well as the free energy. We find that for a surprisingly large 
range of p the results from the FFA agree very well with the results from a reference simulation in the dual 
formulation of the Z3 spin model which is free of the complex action problem. 


To appear in Physics Letters B. 


Introductory comments: 

It is well known that at finite density Monte Carlo simulations of many lattice field theories are plagued by the 
complex action problem (or sign problem). At non-zero chemical potential p the action S acquires an imaginary 
part and the Boltzmann factor e~ s cannot be used as a probability weight. Several different strategies to overcome 
the complex action problem have been proposed over the years, among them the density of states (DoS) method, 
which is the topic of this letter. 

Originally the DoS method was introduced to lattice field theories in [I] and for some recent applications and 
a critical review see, e.g., HE). The key problem of DoS techniques is that the density of states p varies over 
many orders of magnitude. When applying DoS to finite-^ problems the density p is integrated over with a highly 
oscillating factor and the frequency of the oscillation increases exponentially with p. Thus the numerical challenge 
is twofold: The density p has to be determined over many orders of magnitude and this determination has to be 
very precise since p is probed with a highly oscillating factor. 

An interesting new approach to the accuracy problem in the DoS for lattice field theories was presented in 
[4]. The idea is related to a proposal by Wang and Landau |5) and divides the variable of the density p into 
small intervals where restricted Monte Carlo simulations are performed to determine p in that interval. This 
leads to exponential error suppression [4] and was shown to allow for a precise evaluation of observables in pure 
gauge theory and in SU(2) gauge theory with a Polyakov loop source 0]. Furthermore, in [6] the density and the 
complex phase were studied in the Z 3 spin model with chemical potential, which is a simple effective theory for 
the center degrees of freedom of QCD j7]. This model is particularly interesting since it can be mapped to a dual 
representation without complex action problem where precise Monte Carlo simulations at finite p are possible [8]. 

In this letter we develop the DoS for the Z 3 spin model further (functional fit approach; FFA), using again the 
dual simulation |8j as our reference data. We evaluate observables related to the quark number density, the quark 
number susceptibility as well as the free energy and show that the corresponding dual results can be reproduced 
for a surprisingly large range of the chemical potential. First results with the FFA were presented in [9]. 
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Definition of the model and the density of states: 

The Z 3 spin model in an external field with strength k, a chemical potential p/3 (in units of the underlying QCD 
inverse temperature /3) and an effective temperature parameter r is described by the action 


S[P] = £ 

X 


3 

tJ2(p:p*+> 

. V=\ 


2 ) + ne^{P x - 1 ) + Ke~^(P* - 1 ) , 


(1) 


where the dynamical degrees of freedom are the spins P x £ Z 3 = {1, e i27r//3 , e — i27r / 3 }, living on the sites x 
of a 3-dimensional lattice with periodic boundary conditions. The action is normalized such that S[P] = 0 if 
P x = l\/x. The partition function is obtained by summing the Boltzmann factor over all configurations P, i.e., 
Z = J]{p} The first and second derivatives of I 11 Z with respect to /x/3 have the interpretation of the quark 

number density and the quark number susceptibility. It is obvious that for finite chemical potential, p/3 yf 0, the 
action (jT|) has a non-zero imaginary part and the model has a complex action problem. 

We remark at this point that for the Z 3 spin model p does not play the role of a chemical potential in the 
strict sense, since there is no conserved charge it couples to. Still, we refer to p as the ’’chemical potential" due 
to its origin from the proper chemical potential of QCD and due to the fact that it generates the complex action 
problem in exactly the same way as in QCD. The connection of the effective Z 3 model to QCD is obtained from a 
strong coupling expansion (see, e.g., [10]) combined with a hopping expansion of the fermion determinant. This 
connection also provides the physical interpretation of the spins and of the parameters: The spins P x correspond 
to the Polyakov loops of the underlying lattice QCD formulation, restricted to the center values of SU(3), i.e., 
the group Z 3 . The strong coupling expansion identifies r as a parameter which is an increasing function of the 
temperature T of the underlying lattice QCD theory, while the hopping expansion shows that k is proportional to 
the number of quark flavors and a function of the QCD quark mass m which decreases with increasing m. 

For a convenient notation we introduce abbreviations for the total numbers of spins pointing in each of the 
three possible directions as iV 0 [P] = Yh x & {Px> 1) and N±[P] = $ ( P x , e ±i2,r / 3 ). Obviously W 0 [P] + 1V + [P] + 

iV_[P] = V, where V denotes the lattice volume, i.e., the total number of sites of the lattice. Using these we 
split the action into real and imaginary parts, i.e., S[P] = <S[r[P] +iSi[P], with 


3 

— 2) + k3(N 0 [P] — V) cosh(/z/3) , Si[P] = kV% sinh(p/3) AN[P] , (2) 
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where we have defined AiV[P] = iV+[P] — 1V_[P] £ {— V, — V + 1,..., V}. The variables A+[P], iV_[P] and 
A7V[P] behave under complex conjugation of all spin variables P x —> P* as N + [P] — > 1V_[P], ^V_[P] — > iV+[P], 
A7V[P] — ► — AJV[P], and consequently 5[P*] = >S[r[P] — iS/[P] = 5[P]*. Exploring this relation the partition 
sum can be written as 

Z = ^ eS[P1 = eSRlP] eiSl[P] = H eSRlP] cos ( S i [ p D = 5Z eSfl[P1 cos («V / 3sinh( i j0) AA[P]) . (3) 

{P} {P} {P} {P} 


We now define a weighted density of states (5 here denotes a Kronecker-delta), 

P(d) = e SR[p] 6(d — AJV[P]) , d=-V,-V + 1, .... V - 1, V . (4) 

{ p} 


Exploring again the symmetry properties of the action one trivially finds that p{d) is an even function of d, i.e., 
p(—d) = p(d). Using the density of states, the partition sum ([3]) can be written as 

V 

Z = ^ p(d) cos (^kV3 sinh(///3) ctj . (5) 

d=-V 


Expectation values of observables O(AN) which are a function of AN are given by 


<°> = z £ p(d) ^ cos ^t\/3 sinh(/i^) <Hj Oe {d) + fsin ^tv / 3sinh(/z/3) dj Oo{d) 
^ d=-V 
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where Oe and Oo denote the even and odd parts of O(AN). 

The partition sum © and the expectation values ([ 6 ]) are obtained by summing the density p(d) with the 
factors cos («;v / 3sinh(/j,/3) d) and sin (Ki/3sinh(/t/3) d). While the density p(d) is strictly positive, these factor 
are oscillating with d and the frequency of oscillation increases exponentially with the chemical potential /t/3 and 
linearly with the strength parameter n. Thus for larger values of /i/3 (or n) the density p(d) has to be computed 
very accurately. This is how the complex action problem manifests itself in the density of states approach. 

Computing the density of states : 

For the numerical computation we parameterize the density of states p(d) as 


Ml , Ml \ 

P(d) = II e~ a i = exp f - ^2 a i ) > d = ~ v - ~ v + L - v - v > ( 7 ) 

3 =o ' 3=0 7 

with real parameters a 3 . Note that this parameterization is exact in the sense that it contains V + 1 parameters, 
precisely the number of independent degrees of freedom p{d) has (remember that p(d) is an even function). We 
also remark, that an overall normalization of p{d) can be chosen freely, since it cancels in the expectation values 
© Here we choose the normalization p(0) = 1, which corresponds to setting ao = 0. 

For the calculation of the coefficients a 3 we define restricted expectation values {{0)) n ( A), n = 0,1,... V — 1, 
which depend on a free parameter A, 

<(0»n(A) = y29 n (AN[P}) e s «W e XAN ^ 0(AN[P}), Z n (X) = V 0 n (AN[P]) e s ^ e x (8) 
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otherwise 


for n = 1, 2 ,... V — 1 . 
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Varying the parameter A in ([ 8 ]) probes the density in the interval set by 6 n {d) and the response of the system to 
changing A can be used to determine the parameters of p{d) in that interval. In the restricted expectation values 
([ 8 ]) only real and positive weight factors appear, such that they can be evaluated with a restricted Monte Carlo 
strategy which we will discuss below. 

In particular we are here interested in the observable O = AN, and now use the density of states p(d) in the 
form of ([T]) to evaluate the restricted expectation values ((AN)) n (X). A straightforward calculation gives 


« AAr »oW = 


«A3V»„(A) — n = 
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( 10 ) 


The right hand sides are simple functions of A: They are monotonically increasing (the derivatives with respect to 
A are easily shown to be positive) and for n > 1 have a single zero (limA _ > ._ 00 is negative, lim^+oo is positive). 
Examples for different n are shown in Fig. [l] 

Using Monte Carlo simulations we can evaluate ((AN)) n (X)—n for different values Ai, i = 1, 2,... N\ (typically 
N x = 0(10)) and fit the results according to the right hand sides of ( [To] ) . The one-parameter fit for ((AN)) 0 (X) 
determines the first non-trivial coefficient a± (remember that we chose the normalization ao = 0). The fit value 
for ai can then be inserted in the right hand side of (10) for n = 1 such that with another one-parameter fit of 
((AA r ))i(A) — 1 we can determine < 22 , which in turn is then inserted in the fit function for (( AN)) 2 (X ) — 2, which 
gives <23 from a one-parameter fit, et cetera. Using this sequence of fits we can determine all coefficients aj from 
fits of the Monte Carlo data with simple functions that depend on only a single parameter (compare Fig. |T]). Thus 
we refer to our approach sketched in this letter as "functional fit approach" (FFA). 

We expect that this method has smaller statistical errors for the same numerical effort than iteration methods 
such as the LLR [4i[ 6 ] or some general root finding procedure for the ((AN)) n (X)—n (which also provides iterative 
equations for the a 3 ). The advantage of the FFA comes from the fact that all Monte Carlo data, i.e., the results 
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Figure 1: Monte Carlo results (symbols) for ((AN)) n (\)—n for n = 0, 2, 10, 50, 100, 300, 600, 900, 950, 
990, 995 and 998 as a function of A. The data are for a 10 3 lattice with parameters r = 0.16, n = 0.01 and 


IJ./3 


1.0. The full curves are the fits with the functions on the rhs. of (101. The resulting values for the 


corresponding coefficients a 3 are given in the legend. 


for (( AN)) n (Xi ) — n at all values A^ of the parameter A, are used to determine the coefficients aj. Furthermore, 
it can be seen that possible instabilities from the statistical errors of the Monte Carlo data are minimized here. 


Restricted Monte Carlo: 


The FFA variant of the DoS method described here is based on fitting the Monte Carlo data for the restricted 
expectation values ((AN)) „(A) as defined in ([ 8 ]). For this purpose we first need to generate an initial configuration 
P of the spin variables, such that the constraint AJV[P] G {n — l,n, n+ 1} is obeyecQ Such a configuration 
can easily be constructed by hand, but of course needs to be equilibrated before taking measurements. For this 
and the subsequent computation of observables a slightly modified Monte Carlo update can be used. It contains 
an additional restriction which rejects trial configurations that violate the constraint A N[P] G {n — l,n,n + 1}. 
The acceptance rate is very good throughout and only for n very close to the maximum value of n = V (i.e., the 
cases n = V — 2, V — 1, V") we observe a drop in the acceptance rate. In principle it is easy to compute p(d) for 
these largest values of d exactly with a low temperature expansion. However, since for the values of d where the 
quality of the Monte Carlo data decreases p(d) is already very small, we simply use the data as we obtain them 
from the simulation. In this letter we show results for lattice volumes of 10 3 and 16 3 and in both cases used 10 6 
equilibration sweeps and a statistics of 10 6 measurements separated by 100 sweeps for decorrelation, and all errors 
we display are statistical errors. 

In Fig. [l]we show the Monte Carlo results (symbols) for ((AN)) n (A) —n with n = 0, 2, 10, 50, 100, 300, 600, 
900, 950, 990, 995 and 998 for several values of A in the interval [—7, 7], The data were generated on 10 3 lattices 
for the parameter values r = 0.16, k = 0.01 and p/3 = 1.0. The figure demonstrates that the Monte Carlo data 
show the expected simple behavior as a function of A and can easily be fit (we use a standard \ 2 procedure) with 
the functions given in the right hand sides of ( 10 ). In Fig. |Tj the results of the fits are shown as full curves and 


lr This is the constraint for n > 0. The modification to the n = 0 case is trivial and we omit the discussion of this special case. 
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Figure 2: Results for the logarithm of p/d) as a function of d from a 16 3 lattice for different values of p/3. 
The data we show in the lhs. plot are for r = 0.16, k = 0.01. On the rhs. we use r = 0.178 and k = 0.001. 
The error bars are smaller than the line-width. Note the different vertical scales for the two plots. 


obviously describe the numerical data very well. 

Once the coefficients a,j are determined from the fits, we can build up the density of states p/d) as given in 
0 Results for the density p/d) at different values of p/3 are presented in Fig. [2] The data we show in the lhs. 
plot are for 16 3 lattices with r = 0.16, n = 0.01, while in the rhs. plot r = 0.178, k = 0.001 were used. It is 
remarkable that the range of the values for p/d) strongly depends on the parameters, including also p/3. This is 
due to the fact that the weighted density we use here also includes the Boltzmann factor e Sfl . 

Results for physical observables: 

Having determined the density p/d) we can finalize the calculation and evaluate the expectation values of observ¬ 
ables using ([6|. We consider (M — M *) and xm-m* = ((M — M*) 2 ) — ( M — M*) 2 , where 

M = = No[P] - ±(N + [P]+N-[P])+i^-(N+[P]-N-[P]) , (11) 

X 

such that M — M* = i-\/3(lV + [P] — TV_(P]) = iy/3d. Thus for the evaluation of (M — M*) only an odd part 
appears in (|6| and thus this expectation value is real. 

In Fig.Brwe show the results from the FFA for (M — M*) and Xm-m* on a 16 3 lattice as a function of p/3 (red 
circles). We display results for two sets of parameters, r = 0.16, k = 0.01 on the lhs., and r = 0.178, n = 0.001 
(rhs.). The data are for 16 3 lattices and the statistics is the same as used for the density discussed above. As 
reference data in Fig. [3] we also show the results from a dual simulation. For the r = 0.178, k = 0.001 data (rhs.) 
the FFA results agree well with the dual results for all values of p/3 we show. However, for r = 0.16, k = 0.01 
(lhs.) the dual and the plain FFA data disagree for p/3 larger than 2. This discrepancy can be attributed to the 
fact, that for this parameter set the sign problem is much harder than for r = 0.178, k = 0.001, as can, e.g., be 

seen in plots for the expectation value of the phase of the action (Fig. 1 of [ITj). 

The numerical problems can be related to small statistical fluctuations of p{d) around its exact values [6|. 
However, it is straightforward to smoothen these local fluctuations by using a fit of p/d) and then computing the 
observables ([6]) with the fitted p/d). For the fit we use a polynomial in d 2 (note that p/d) is an even function and 
lnp(0) = 0) for the logarithm of p/d), i.e., In p/d) = Yln =l c nd 2n . In Fig. [ 3 ] we also display the results obtained 
from the fitted p/d) with TV = 15 (blue triangles) and find that we obtain a much larger range in p/3 where FFA 
and dual simulation agree also for the r = 0.16, n = 0.01 data. 

We stress that using a fit to smoothen p/d) is of course not a fundamental ingredient of the FFA, and increasing 

the statistics when determining p/d) will also improve the results. However, the source of the error is very clear: 

For large p/3 the density p/d) is probed by the rapidly oscillating factors in ([ 5 ]) and ([6]) and the small fluctuations 


5 












Figure 3: Results for the physical observables (M — M*) (top row of plots) and Xm-m* (bottom) on 16 3 
lattices as a function of /fob We use two sets of parameters, r = 0.16, k = 0.01 on the lhs., and r = 0.178 
and k = 0.001 on the rhs. We show results from the FFA algorithm, the FFA algorithm combined with a fit 
of p{d) and for comparison also the results from a simulation in the dual representation. 


in p{d) become dominant. On the other hand p{d) is a smooth function, such that a fit is a much more cost 
efficient method than a drastic increase of the statistics. 

The influence of small fluctuations of the density p(d) is studied in Fig. [4] In this plot we zoom into the lower 
values of d (note that d runs from 0 to V = 16 3 = 4096), and show as a function of d the density p(d), the os¬ 
cillating factor sin(K\/3sinh(/z/3)cJ) and the cumulative sum S(d) = —2y/3/VZ'^2j_ 1 p(j)sin(Ky/Zsinh(fj,^)j)j, 
which for d = V sums up to (M — M*)/V. Thus studying the cumulative sum S(d) as a function of d shows 
how (M — M*)/V is built up. 

For the smaller value pfd = 1.0 shown in the lhs. plot, we find that S(d) shows no sizable fluctuations and 
approaches its asymptotic value without major fluctuations. For p/3 = 3.0, however, we observe that S(d) strongly 
picks up the now faster oscillations from sin(re\/3sinli (p/3)d) and crosses the value S(d) several times. This implies 
that large cancellations are necessary to reach the asymptotic value and small fluctuations of the density p(d) have 
a large impact. These fluctuations are suppressed when fitting the density as discussed above. 

As our final observable we consider the free energy F(pf3) defined as F{pj3) = \n.(Z{pj3)/Z(Q)). In Fig. [5] we 
show the results for the free energy at two values of the couplings and again compare the FFA method with data 
obtained from a dual simulation. As for the other observables we find excellent agreement of the FFA results with 
the data from the reference simulation in the dual approach. It is remarkable that for the free energy this excellent 
agreement is achieved already without fitting the density, which was necessary for the bulk observables. 
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Figure 4: Buildup of the contributions to (M — M*)/V for V = 16 3 , r = 0.16 and k = 0.01 at /i/3 = 1.0 
(lhs.) and /i/3 = 3.0 (rhs.). See the text for a discussion of the plots. 




Figure 5: Results for the the free energy F(p,/ 3) from the FFA method (circles) and from the dual simulation 
(crosses). The data are for V = 10 3 with r = 0.16, n = 0.01 (lhs.), and r = 0.178, k = 0.001 (rhs.). 

Summary: 

The key issue in the application of the density of states method is to determine the density p(d) as precise as 
possible, since in the expectation values ([6]) p(d) is summed over with a highly oscillating function. Thus it is 
necessary to optimize the strategy for the computation of p(d) in every possible way. In these notes we test a 
strategy (the functional fit approach (FFA)) where in restricted Monte Carlo simulations on small intervals on d 
the response of the system to a free parameter in the Boltzmann factor is evaluated. The response is given by a 
known function which we fit to the Monte Carlo data to determine the parameters of the density p{d). 

The results from the DoS calculation for (M — M*), Xm-m* and the free energy F(p(3) are compared to 
reference data obtained in a dual simulation. We show that when using a fit for the density of states p(d) the 
observables agree very well with the dual results on a surprisingly large range of /i/3 values up to /i/3 « 4. 
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